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PREFACE 


This survey paper was presented to the AGARD Fluid Dynamics Panel in October 
1975 to aid the panel in formulating plans for a future conference on theoretical and 
computational aerodynamics. It describes in part programs in theoretical and com- 
putational aerodynamics in the United States. It details only those aspects of those 
programs that relate to aeronautics and not the related development effort involving 
application of aerodynamic theory to such problems as planetary entry and atmospheric 
and oceanographic modeling. The role of analysis at various levels of sophistication is 
discussed as well as the inverse solution techniques that are of primary importance in 
design methodology. The research discussed is divided into the broad categories of appli- 
cation for boundary layer flow, Navier-Stokes turbulence modeling, internal flows, two- 
dimensional configurations, subsonic and supersonic aircraft, transonic aircraft, and the 
Space Shuttle. A survey of representative work in each area is presented, with the ex- 
ception of turbulence modeling techniques, which were covered in a recent AGARD 
presentation. 
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INTRODUCTION 


The role of computational aerodynamics is to aid 
in advancing technology of long-range aircraft sys- 
tems, helicopters and rotorcraft, general aviation air- 
craft, and maneuverable fighter aircraft and missiles 
with the generation of design and analysis methods 
that provide improved definition and lower costs in 
developing such vehicles. The development of new 
aerodynamic concepts evolves from a clear under- 
standing of the physics of flow fields. To obtain this 
understanding with the present state of the art, an 
appropriate balance between theory and experiment 
must be maintained at all research levels. Although 
new design and analysis methodology may take the 
form of computational techniques, its development 
also rests squarely on the dual foundations of theory 
and experiment. 

Fluid flow research has generally emphasized the 
complementary roles of theory and experiment. For a 
number of years, the experimentalists directed most 
of their effort toward determining aerodynamic char- 
acteristics for point design and obtaining systematic 
data for configurations and flow phenomena not yet 
amenable to analysis. Experimental programs for vali- 
dating a theoretical method or for providing data for 
the generation of a theoretical model were few. Theo- 
retical methods were based primarily on linear, in- 
compressible, inviscid flow equations or on the 
Prandtl boundary layer equations and were applicable 
to only the simplest configurations. With the formula- 
tion in the late 1940s and 1950s of a number of 
methods for solving approximate linear and exact 
forms of the governing equations for the flows over 
practical geometries, an increased emphasis on experi- 
mental research for theory validation took place. This 
trend has accelerated in recent years, stimulated by a 
rapidly improving analysis capability, which may be 
attributed to growth in the speed and storage capac- 
ity of electronic computers as well as the develop- 
ment of improved numerical solution techniques. 

In 1946 John von Neumann described wind tunnel 
testing as follows (ref. 1): 

Indeed, to a great extent, experimentation in fluid 

dynamics is carried out under conditions where the 


underlying physical principles are not in doubt, where 
the quantities to be observed are completely deter- 
mined by known equations. The purpose of experi- 
ment is not to verify a theory but to replace a com- 
putation from an unquestioned theory by direct 
measurements. Thus wind tunnels, for example, are 
used at present, at least in part, as computing devices 
to integrate the partial differential equations of fluid 
dynamics. 

As computer simulations become more feasible be- 
cause of rapidly improving computer capability, the 
need for as much wind tunnel testing will decrease. 
The dual partnership of theory and experiment con- 
tinues; however, their roles in relation to aeronautics 
programs are beginning to change with the expanding 
role of the computer that has occurred because of its 
increasingly accurate simulations, as indicated in 
figure 1 . The requirement for point design experimen- 
tal data has steadily diminished as numerical methods 
have improved. Limitations imposed by most experi- 
mental programs, where only a few configurations 
can be tested over a limited range of conditions, are 
being circumvented to an increasing extent by empha- 
sizing tests specifically designed to check theory at a 
few critical points or to generate and improve theo- 
retical flow models. Once an analysis is proven ade- 
quate, it provides an interpolative and extrapolative 
capability for handling other geometries or flow con- 
ditions for which test data do not exist. As the cost 
of wind tunnel testing increases, computer simula- 
tions are expected to reduce the time and cost of 



Figure 1 .-Comparison of accuracy of wind tunnel testing 
and computer simulations. 
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developing new prototype vehicles. In summary, 
then, computational fluid dynamics has the following 
advantages over wind tunnel testing: 

(1) It provides the capability of simulating flow 
impossible or impractical to simulate in a ground- 
based experimental facility. 

(2) It provides more accurate simulations of flight 
aerodynamics than wind tunnels can provide. 

(3) It reduces the time and cost of obtaining flow 
simulations necessary in the design and analysis of 
new prototype aerospace vehicles. 

In the 1980s computer methods and wind tunnel test- 
ing will perform complementary functions. 

New numerical techniques and large computers 


have also allowed solutions of problems that, until a 
few years ago, could not have been attempted even 
with the most gross assumptions. Solutions of the 
Navier-Stokes equations, or slightly simplified ver- 
sions, for viscous flow over multielement airfoils, 
shock/boundary-layer interactions, various jet-ex- 
haust and free shear flows, and for a variety of two- 
dimensional geometries where flow separation is an 
important feature, indicate a new level of technology 
that will mature swiftly with the development and 
use of computer systems such as the CDC STAR or 
the ILLIAC IV. With the use of such or even more 
advanced computers, a similar trend may take place 
in the application of numerical simulation to flows 
over complex aerospace configurations. 



METHODS AND STATUS 


The basic process used in theoretical and computa- 
tional aerodynamics is the solution of a set of govern- 
ing partial differential equations, which are known, 
under constraints imposed by the initial and bound- 
ary conditions specified as appropriate to the prob- 
lem at hand. The representations (mathematical 
model) of these equations and boundary conditions 
of the physical problem are determined by the choice 
of numerical method, computer language, and by the 
machine itself. A computer code is then developed to 
provide a solution to the modeled problem. Results 
of calculations from the code are compared with ex- 
perimental data to determine the validity of the 
representation. Experimental data are also needed to 
provide empirical inputs for the computer code, and, 
as computations are becoming available for more 
sophisticated flow fields, researchers are able to deter- 
mine from computed solutions areas in which experi- 
ments will be particularly helpful in understanding 
aerodynamic phenomena. Thus the cycle of experi- 
ment and analysis illustrated in figure 2 continually 
advances both efforts. 

The time-dependent Navier-Stokes and energy 
equations are generally accepted as a valid set of gov- 
erning equations for aerodynamic flows without 
chemical reactions. In 1972 Bradshaw stated (ref. 2), 

In turbulence studies we are fortunate in having a 
complete set of equations, the Navier-Stokes equa- 
tions, whose ability to describe the motion of air at 
temperatures and pressures near atmospheric is not 
seriously in doubt. (It is easy to show that the smallest 
significant eddies are many times larger than a molec- 
ular mean free path.) We are unfortunate because 
numerical solution of the full time-dependent equa- 
tions for turbulent flow is not practicable with present 
computers. 

Because the Navier-Stokes equations are highly 
nonlinear, their exact solution is virtually impossible. 
Useful solutions in aerodynamic theory have been 
made possible by the application of simplifying 
approximations of one form or another to the equa- 
tions and their boundary conditions. Some standard 
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Figure 2. -Elements of a computational fluid dynamics 
program. 

approximations, using the * momentum equation 
alone for illustration, follow, in increasing order of 
approximation. The x momentum equation is 


du du 

l ‘\W tu i to, 


dp d / du\ 
dx + d Xj\ dxj 



In the Reynolds model, velocity is assumed to be the 
sum of an average and a fluctuating quantity; i.e., 


therefore, 


P 


du du\ 
dt + Uj ~dXjj 


dp d / du \ 

dx + dxj y ^ dXj J 


2 d / duA d /_dwA 
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In boundary layer modeling, it is assumed, in addition 
to the Reynolds assumptions, that 

u » v 
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and 


algebraic manipulation, this equation becomes 


du By 
0 y ^ dx 


therefore, 


0 - M l)- + ~ + - = MUy - n" 

' - } dx dy dz °° {7 ’ Udx 


In the in viscid linearized model, the right-hand side of 
the preceding equation is neglected; therefore, 


bu __ bu 

pu ^ + pl, a7 


bp 0 / _ du 

bx by V by 



The inviscid nonlinear model assumes p = 0 in the 
Navier-Stokes equations, thus 



+ P 



+ 




bp 

bx 


After definition of a perturbation velocity and some 




Each form of approximation or subdivision 
thereof leads to its own solution techniques with 
comparative advantages and disadvantages with re- 
spect to generality of application, accuracy of results, 
and cost of solution. Early computers were capable of 
computations simulating only the simplest flows. As 
computer storage and speed have increased, the level 
of approximation that could be simulated has in- 
creased. Table 1 lists the principal limitations in appli- 
cation, the status of available computer codes, and 
the pacing items in code development for various 


TABLE 1 .—Stages of Approximation to Governing Equations of Fluid Dynamics 


Characteristics 

Approximations 

Principal limitations 

Status 

Pacing item 

Inviscid linearized 

Viscous and nonlinear 
inviscid terms 
neglected 

Slender configurations; 
small angle of at- 
tack; perfect gas; 
no transonic flow; 
no hypersonic 
flow; no flow sep- 
aration 

2 -dimensional flows in 
1930s; 3-dimen- 
sional flows in 
1960s; used in cur- 
rent aircraft design 


Inviscid nonlinear 

Viscous terms 
neglected 

No flow separation 

2 -dimensional flows in 
197 l ; 3-dimensional 
flows in 1975; early 
stages of application 
to aircraft design 

Code development 

Viscous boundary 
layer 

St ream wise velocity 
gradients ne- 
glected; turbu- 
lence terms mod- 
eled by eddy 
viscosity 

No flow separation; no 
large crossflow; no 
large pressure 
gradients 

2-dimensional flows in 
1971 ; 3-dimensional 
Hows in 1 974; early 
stages of application 
to aircraft design 

Code development 

Viscous time-averaged 

No terms neglected; 

Accuracy of turbu- 

2-dimensional flows in 

Development of un- 

Navier-Stokes 

turbulent momen- 
tum, energy, and 
heat transport 
terms modeled 

lence model 

1975; 3-dimensional 
flows possibly by 
1978; early stages 
of development with 
simple turbulence 
models for 2-dimen- 
sional flows 

proved turbulence 
models 

Complete viscous time- 
dependent Navier- 
Stokes 

Subgrid scale motion 
modeled 

Accuracy of Navier- 
Stokes equations 

Mid-1980s 

Development of ad- 
vanced computer 
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TABLE 2 —Methods of Solution 


Method 

Advantages 

Disadvantages 

Comments 

Explicit finite differences (for 

Easy to code; incorporation 

Large number of nodes and 

Most frequently used in the 

example, MacCormack, 
hopscotch, and DuFort- 
Frankel) 

of turbulence modeling is 
straightforward 

long computation times 

past because of simplicity 

Implicit finite difference (for 

Potential for shorter com- 

Difficult to code; requires 

Advent of large computers 

example, alternating 
direction implicit (ADI)) 

putation time 

matrix reduction; re- 
quires excessive storage 

has made this feasible 

Finite element 

Same degree of accuracy re- 
quires fewer nodes; 
easily adapted to irregular 
geometry; derivative 
boundary conditions in- 
corporated directly 

Data handling complex and 
requires excessive time 
and storage 

Widely used in structures; 
now being investigated 
for fluid problems 

Cubic spline integration , 

Same advantages as finite 
element; handles shock 
waves easily 

Requires tridiagonal matrix 
reduction 

Recently developed tech- 
nique that still has to 
prove itself in application 


stages of approximation to the Navier-Stokes equa- 
tions. 

As noted previously, the equations represent the 
physical problem. The numerical procedure chosen to 
model the equations provides the mathematical repre- 
sentation of this problem in the computer. Listed in 
table 2 are some of the procedures that are currently 
used, and the principal advantages and disadvantages 
of each method. Because the computation cost and 
storage requirement are a direct function of the com- 
plexity of the equation simulated, the simplest 
approximation of the governing equations that is 
applicable to the physical problem should be used. 
Once the chosen equations are cast in a finite differ- 
ence form by using a numerical technique, the solu- 
tion can proceed either by a time marching or a re- 
laxation method. At present, there is no one optimal 
numerical technique. The choice of the numerical 
method must be based on the flow characteristics of a 
particular application and the operating character- 
istics of the available computational facility. In the 
case of numerical simulation of aircraft configura- 
tions— the shuttle configuration shown in figure 3, for 
example -highly integrated and specialized applica- 
tions of several solution techniques are required de- 
pending on the level of detail necessary and the loca- 
tion on the vehicle where the flow is of interest. 



Figure 3. -The combination of numerical procedures used in 
the analysis of a single configuration. Region I: subsonic 
How; must march in time. Region II: supersonic flow; no 
plane normal to body axis with all supersonic flow; use 
space marching mapped coordinates. Region III: super- 
sonic How; planes normal to body axis with all supersonic 
flow do exist; use parabolic Navier-Stokes equations with 
space marching body coordinates. Region IV: subsonic 
How; embedded regions of subsonic flow around wing and 
tail leading edges; use full Navier-Stokes time marching or 
couple potential equations with boundary layer equa- 
tions. Region V: internal How; reacting chemistry must be 
included. Region VI: nozzle exhaust; three-dimensional 
mixing; shock interacting. 


SURVEY 


A survey of representative U.S. work in computa- 
tional aerodynamics is presented in tables 3 to 8 in 
terms of the categories of boundary layer flow, in- 
ternal flows, two-dimensional configurations, sub- 
sonic and supersonic aircraft, transonic aircraft, and 
space shuttle. Highlights in each category are dis- 
cussed and illustrated with specific examples in the 
following sections. The large amount of effort 
directed toward Navier-Stokes turbulence modeling is 
not included in tables in this document because it is 
the subject of a separate AGARD review (ref. 3). 
Some discussion of this effort is presented, however, 
because of its importance to the future of the entire 
computational aerodynamics effort. 

VISCOUS FLOWS 

Viscous flow aerodynamics is discussed prior to 
and apart from the internal flow and configuration 
categories because of its importance and across-the- 
board applicability. Work in this category divides it- 
self naturally into the areas of boundary layer flows 
and Navier-Stokes solutions. 

Boundary Layers 

Numerical solutions for two-dimensional and 
three-dimensional laminar boundary layers over 
simple shapes without large adverse pressure gradients 
can presently be achieved to any desired accuracy. 
With existing models for the Reynolds stresses, turbu- 
lent boundary layer flows can be calculated with rea- 
sonable accuracy. There has also been some success in 
computing 2-dimensional flow through small em- 
bedded separation regions using the boundary layer 
equations; however, 3 -dimensional flow near separa- 
tion is not fully understood. Although prediction 
methods for laminar and turbulent flow are accurate, 
little is known about the prediction of transition and 
relaminarization. Current research is concentrating in 
these areas and in developing more efficient proce- 
dures. A sampling of this research is given in table 3. 
No distinction is made in the table between codes 


that are for research only and those that are or will be 
available as production codes. 

There have been several survey papers on solution 
procedures for the boundary layer equations. Also, 
books are available describing some of the methods in 
detail. Second-order accurate, implicit and explicit 
finite difference schemes exist for solving steady 
two-dimensional and axisymmetric perfect gas flow. 
The method of weighted residuals, the matrix integral 
method, and higher order finite difference methods 
have also been applied with success. Optimal spacing 
of computational nodes normal to the wall boundary 
plays a vital part in retaining the desired accuracy- 
while minimizing computer time for turbulent flows. 
As more efficient techniques are developed, they can 
be used in unsteady, three-dimensional, and real-gas 
flows. Some two-dimensional techniques have been 
extended to the unsteady flow case. Additional work 
on the appropriate transformation to use with the 
unsteady boundary layer equations is needed. 

For steady three-dimensional flows, difference 
schemes have been developed and problems have been 
solved for various initial and boundary conditions. 
One problem appears to be the development of a 
general code for solving a variety of three-dimensional 
flows. Such a development requires further evaluation 
of coordinate systems, transformations, and more 
accurate and better ways of handling the required 
inviscid flow data. Flexibility needs to be added into 
the codes so that the various difference schemes can 
be automatically used as required to satisfy zones of 
dependence. There appear to be no new problems in 
extending three-dimensional numerical techniques to 
unsteady flows, other than computer storage require- 
ments, unless reverse flow is encountered. Some 
three-dimensional solution procedures have been ex- 
tended to turbulent flows with no more difficulty 
than encountered for two-dimensional extension to 
turbulent flows. One implicit code used by a number 
of investigators employs a second-order differencing 
scheme that is conditionally stable for reverse cross- 
flow. It has an option for either two-layer eddy vis- 
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TABLE 3 —Boundary Layer 


Investigator 

Method 

Application 

Code status 

J. Adams, 

Implicit finite difference 

; Yawed airfoil 

Operational for 3 dimensions 

H. Dwyer 

Finite difference, hybrid 
scheme 

, Spinning cone 

Operational for 3 dimensions 

J. Nash 

Differential method 

! Airfoil sections 

Operational for 3 dimensions 

A. J. Baker 

Finite element 

Combustors; duct Hows 

Operational for 3 dimensions 

K.C. Wang 

Mixed one- and two-step 
finite difference 

Blunt body 

Operational for 3 dimensions 

A. Wortman 

Differential operator tech- 
nique 

Blunt body; swept and 
tapered wings 

Operational for 3 dimensions 

R. T. Davis 

Implicit finite difference 

Sharp and blunted cones 

Operational for 3 dimensions 

F. Blottner 

Finite difference 

Analytic bodies 

| Operational; plans include 
incorporation of flexible 
integration direction 

R. M. Kendall and M. Abbett 

Implicit solution: finite dif- 
ference and splines 

Windward surface with rea* 
gas effects 

Operational; plans include 
coupling with 3-dimen- 
sional inviscid code 

T. Cebeci 

Keller “box” method 

Perfect gas analysis for wings 
of arbitrary shape 

Operational with eddy vis- 
cosity model; plans in- 
clude coupling with 3- 
dimcnsional inviscid 
code 

J. Harris 

Implicit finite difference 

Analytic bodies; option for 
viscosity model 

Analysis operational for 3 
dimensions; plans include 
extension to general 
geometry 

A. J. Baker 

Finite element 

Reentry vehicle in real gas 

Operational for 2-diinen- 
sional supersonic How; 
plans include extension 
to subsonic flow 

M. Frieders and C. Lewis 

Implicit finite difference 

Sharp and blunt cones with 
real gas effects 

Operational for 3 dimensions 


cosity or mixing length turbulence models. Investiga- 
tions have been made concerning reduction in nodal 
spacing to optimize computer time and minimize 
storage. Some typical results presenting calculations 
using the two different turbulence models and two 
nodal spacings are shown in figure 4 (from ref. 4). A 
typical case of 101 points in the normal direction 
requires 300 s and 70 000 octal storage locations on a 
CDC 6600 machine. An example of the capability of 
the current analysis at hypersonic speeds to handle 
real-gas effects on a blunted cone at angle of attack is 
shown in figure 5 (from ref. 5). The good agreement 
on the leeward side is not surprising because of the 
relative small angle of attack of the surface. 

Navier-Stokes and Turbulence Modeling 

For viscous flows where the boundary layer 
assumptions are invalid, the full time-dependent 
Navier-Stokes equations must be solved. Such solu- 


tions are presently obtainable for only a restricted 
class of these flows because of limitations in the finite 
difference computations. Turbulent flow is particu- 
larly difficult to calculate. At a point in space, turbu- 
lent flow is irregular and random. It is a three- 
dimensional phenomenon composed of different 
structures having a wide range of time and length 
scales. To compute these irregularities, accurate time- 
dependent computational techniques must be used. 
Numerical stability and accuracy restraints on the 
allowable time step are imposed by the small-scale 
dimensions required to prevent errors in the instan- 
taneous flow field from growing. These requirements 
for computer speed and storage capacity are beyond 
available and projected computer capabilities; thus it 
is necessary to resort to turbulence modeling. 

The first step in turbulence modeling is the elimi- 
nation of the small-scale structure of the dependent 
variable by an averaging process (ensemble, space, or 


SURVEY OF VISCOUS FLOWS 


9 



(a) 



Experimental data: 
o u/Q e 
□ v/Q e 
a TfT e 


Computer data: 

x Two -layer eddy viscosity 
turbulence model 

Mixing length turbulence model, 

301 data points 

__ Mixing length turbulence model, 
21 data points 


Figure 4. -Experimental data and computer calculations for 
two different turbulence models at two nodal spacings 
(from Harris and Morris, ref. 4). M ^ = 1.8; Re^ L = 25 X 
10 6 ; 7 = 1.4; Q e = (u e 2 + v e 2 ) 1/2 - (o) <t> = 0°.' (ft) <l> = 
135°. 



Storage and processing on the CDC 6600 


Number of 

Central processing 

Storage 

Test case 

data points 

unit mesh point, s 

time, s 

301 

0.0123 

170 000 

900 

21 

.0123 

60 000 

60 


time). The new averaged dependent variables can be 
resolved with finite difference methods using mesh 
spacing and time intervals compatible with present 
computer capability. The nonlinear terms in the basic 
equations introduce averages of products of the 
small-scale portions of the dependent variables that 
must be related to the averaged values (the “closure” 
problem). This “closure” is the second step in the 
modeling process. Here much emphasis is placed on 
comparisons with experimental data because the 
modeling is empirical. 

One method under investigation, suitable for 
steady-state mean flow fields, is called statistical 
theory of inhomogeneous turbulence, where averag- 
ing is performed at a point in space over a long period 
of time compared with the time scales of the largest 
turbulent eddies. The closure is achieved by eddy vis- 
cosity modeling. An example of the effectiveness of 
the eddy viscosity turbulence model in transonic flow 
as applied to a circular arc airfoil is given in figure 6 


(from ref. 7). The model is shown to be quite effec- 
tive when compared with the inviscid calculation in 
the region of attached flow or when the separation is 
modest as occurs in trailing edge separation. With 
shock-induced separation it is clear that the model 
used is inadequate although the major features of the 
flow are qualitatively predicted. Figure 7 (from ref. 
8) compares the experimental results of a well- 
documented experiment in which a shock wave was 
used to separate a turbulent boundary layer at mach 
7.2. The predicted quantities were calculated with a 
statistical turbulence model of the Navier-Stokes 
equations. Although the general features of the data 
are computed, there are significant differences that 
require improvement. The dashed lines represent an 
improved model in which the data were analyzed to 
change the principal modeling constants. Some im- 
provements occur, especially with regard to the up- 
stream influence, but the algebraic models are most 
limited and much needs to be done to improve the 
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O Experimental data 

— Viscous calculation 

— Numerical inviscid calculation 




(x - x 0 )/5 0 


O Experimental data 

— Baseline model 

— Modified model 


Figure 6. -Numerical solution of circular arc airfoil using 
Navier-Stokes equations with eddy viscosity turbulence 
model (from Deiwert, ref. 7). (c^* = critical pressure co- 
efficient.) (a) Trailing edge separation. M ^ = 0.775; 

oo = 2 X 10 6 . (b) Shock-induced separation. = 
0-7&6; Re c , = IQ 7 . 


Figure 7. -Pressure, skin friction, and heating calculations 
for oblique shock wave/turbulent boundary layer inter- 
action at mach 7.2 (from Marvin et al., ref. 8). (< Cp - 
skin friction coefficient; = heat transfer coefficient.) 
(a) a = 7.5°. 
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□ Experimental data 
— Baseline model 

__ Modified model ^ 

(b) 

Figure 7 (concluded).- (b) a = 15°; S = separation;/? = at- 
tachment. 

modeling for this type of a complex flow field. 
Another example of the effectiveness of the eddy vis- 
cosity model is shown in figure 8 (from ref. 9) for the 
case of a normal shock wave turbulent boundary 
layer interaction at mach 1.5. The numerics ade- 
quately predict the skin friction up to separation but 
underestimate skin friction downstream of reattach- 
ment. However, the length of the separated flow 
region is fairly closely predicted. Also, although the 
simulation considerably underpredicts the height of 
the reverse flow region, it does provide a good quali- 
tative description of the flow by identifying the ex- 
tent of the shear layer and the location of maximum 



x/8 u 



o Profile data 

— Computation 

— Preston tube 



o Experiment 
□ Computation 



(b) 

Figure 8. -Normal shock wave/turbulent boundary layer in- 
teraction in a circular channel at mach 1.5 (from 
Mateer, Brosh, and Viegas, ref. 9). (tf) Comparison of 
wall pressure and skin friction. ( b ) Normal shear stress 
distribution; x/s = 4. 
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shear. Ultimately, a prediction method may have to 
drop the eddy viscosity concept and calculate the 
Reynolds stresses directly. 

An alternative method, sometimes called “turbu- 
lent simulation,” is based on averaging over space 
volumes that are smaller than the largest turbulent 
eddies but much larger than the smallest eddies. The 
volume averages retain their time dependence but 
possess scales resolvable by current computer tech- 
niques. Turbulent simulation is expected to become 
the basis for practical turbulence computation be- 
cause inherent in its formulation is its direct evalua- 
tion of the large eddies that are characteristic of a 
particular flow and its modeling of the small “sub- 
grid” scales that are thought to be universal in 
character (fig. 9). 

INTERNAL FLOWS 

Aircraft with airbreathing propulsion require a 
high degree of engine/airframe integration to achieve 
optimized performance. Although the complex flow 
interactions resulting from integrated aircraft configu- 
rations present a formidable challenge to analysis, the 
effects of these interactions must be determined to 
produce an effective aircraft design. Representative 
research in this area is listed in table 4. 



Figure 9. -Goal of turbulence modeling with projected 
computer capability. Subgrid-scale turbulence modeling is 
used to close the system of equations. 


Inlets 

Inlet flows represent one of the more difficult 
challenges in the development of computational 
techniques. For the most part, inlets have been 
developed through a combination of empiricism and 
many hours of wind tunnel testing. Analytical tech- 
niques are needed for the detailed design of internal 
contours, mass removal systems, and for future 
aircraft to blend the propulsion system into the 


TABLE 4 -Internal Flows 


Investigator 

Method 

Application 

Code status 

D. A. Oliver and P. Sparis 

Finite difference, inviscid 

! Turbomachinery 

Operational 

J. Erdos and E. Alzner 

Finite difference, inviscid 

Turbomachinery 

Operational 

L. L. Presley 

Finite difference, shock- 
capturing, inviscid 

Supersonic inlet 

Operational 

P. Kutler and V. Shankar 

Finite difference, inviscid 

Integrated ram/scram jet 
nacelles 

Operational 

S. M. Dash and P. D. Del 
Guidice 

Finite difference 

3-dimensional nozzle ex- 
haust flow hydrogen-air 
chemistry 

Operational 

C. du P. Donaldson 

Finite difference, second- 
order closure 

Combustion in scram jets at 
low mach numbers 

Operational 

R. S. Hirsh 

Finite difference (ADI), x- 
marching, viscous 

3-dimensional nozzle exhaust 

Operational for supersonic 
flow; plans include ex- 
tension to subsonic 

H. McDonald and W. R. 
Briley 

Finite difference (ADl),x- 
marching, viscous 

Internal flows 

Operational 

W. S. Llewellen and R. D. , 

Sullivan 

Finite difference, x- 
marching, viscous 

Subsonic wake flow now 
being applied to vortex 
mixing 

Operational 

M. D. Salas 

Finite difference (Mac- 
Cormaek) and floating 
shock fitting, inviscid 

Supersonic inlet 

Operational 
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overall configuration with minimum penalty to, or 
perhaps enhancement of, the basic aircraft perform- 
ance. 

Most of the analytical work in the past has been 
applied to the development of computer codes for 
planar and axisymmetric supersonic flows using the 
method of characteristics. Pressure distributions from 
these codes are used as input to viscous codes 
incorporating engineering models of shock-wave 
boundary-layer interactions to predict the extent and 
magnitude of boundary layer removal. 

Solution of subsonic and transonic internal flows 
is not as well developed as supersonic flows. Ap- 
proaches to subsonic flows have usually resulted from 
incompressible potential flow solutions, with com- 
pressibility corrections being applied. Transonic inlet 
flows have been approached using small disturbance 
theory as well as time-dependent solutions of the 
equations of motion. Significant progress has been 
made in analyzing subsonic duct flows in which 
viscous effects and arbitrary cross section are in- 
cluded in the analysis. Recent developments have 
been along several lines: 

(1) Developing techniques for full three-dimen- 
sional supersonic flows, including corners 

(2) Developing three-dimensional techniques for 
transonic flows 

(3) Developing improved viscous techniques for 
analyzing flows with large adverse pressure gradients 
and mass removal 

(4) Developing analytical techniques that simul- 
taneously analyze the inviscid and viscous portions of 
the flow and include coupling effects 

One method is illustrated in figure 10 (from ref. 10) 
where calculations of pressure ratio in the inlet are 
compared with a solution by the method of charac- 
teristics. 

It should be emphasized that the primary emphasis 
of future efforts should be in developing analytical 
techniques that can be applied to general three- 
dimensional configurations and that viscous effects 
should be included. However, throughout the devel- 
opment of new computational techniques, continual 
comparison with experimental standards must be 
made to verify the accuracy and physical modeling of 
the new techniques. 

Turbo machinery 

The flow field within advanced axial turbomachin- 
ery, e.g., high bypass ratio turbofan engines, is 



Z K 



External flow 



— Method of characteristics 
O Shock-capturing technique 

Figure 10. -Comparison of numerical techniques (from Pres- 
ley, ref. 10). M oo = 3.5; a = 0°. Dashed lines indicate 
region of boundary layer removal. 

characterized as highly three-dimensional and un- 
steady, with transonic to supersonic mach numbers 
and multiple shocks. Conventional analytical methods 
place reliance upon two-dimensional potential flow 
and linearized three-dimensional flow models, and, 
therefore, constrain to some extent the design and 
analysis of new systems. An example of the develop- 
ment of more advanced methods of aerodynamic 
analysis based on numerical solution of the complete 
equations of motion for an unsteady inviscid com- 
pressible flow is the joint effort undertaken by NASA 
Lewis Research Center and by Advanced Technology 
Laboratories, Inc. The programs employ a shock- 
capturing finite-difference method at the interior grid 
points and a reference plane method-of-characteristics 
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procedure at the boundary points. Subsonic axial 
velocity is assumed at both the inlet and discharge 
stations, and either choked or unchoked operation 
can be considered. One program pertains to analysis 
of meridional hub-to-shroud stream surfaces and 
another to axisymmetric blade-to-blade stream sur- 
faces (fig. 11, from ref/ 11). Both programs can 




analyze either one or two blade rows (i.e., rotor 
and/or stator), and the meridional program can 
consider multiple stages. Finite blade thickness, cam- 
ber, and lean angles are included. In the meridional 
program, a steady, axisymmetric solution with cross- 
flow (i.e., swirl) is obtained. In the blade-to-blade 
program, a steady two-dimensional solution is ob- 



(c) 



Experiment {ref. 12) 

-O Present numerical results 
(10 X 8 mesh) 


Figure 11. -Rotor pressure distributions (from Erdos et al., ref. 11). (a) Diagram. ( b ) Experimental 
isobars (ref. 12). (c) Numerical isobars. ( d ) Suction surface. 
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Experiment (ref. 12} 

-o- Present numerical results 
(10X8 mesh) 



Figure 11 (concluded). -(e) Midchannel line, (f) Compression surface. 


tained for a single blade row, and a periodic unsteady 
flow is obtained with a rotor and stator. Numerical 
results to date consist of certain test cases for a single 
transonic rotor blade row (for which experimental 
data are available) and a hypothetical rotor/stator 
combination. The computed results agree favorably 
with experimental data, as shown in figure 11, in 
which the free stream mach number relative to the 
rotor is 1.5. The present computer program produces 
a reasonably good rendition of the experimental data 
after 200 time steps, which required only 140 s of 
CDC 6600 computer time. 

Combustor 

The complex flow field resulting when two or 
more ducts merge, as, for example, the flow produced 
by fuel injectors in a nozzle or the merging flow of a 
series of scramjet nodules, produces a number of 
interacting shocks. A finite difference computer code 
has been developed to analyze the flow of a perfect 
gas inside two-dimensional, axisymmetric, or line- 
source-type ducts. For the axisymmetric case, provi- 
sions to calculate flows with swirl are included. The 
geometries of the ducts are arbitrary, provided that 
the flow remains supersonic. The formation of shock 


waves is automatically predicted, and all shock waves 
and vortex sheets are treated (using the floating shock 
fitting technique) as discrete discontinuities that 
satisfy the appropriate jump conditions. There are no 
restrictions on the possible number of shock waves 
and vortex sheets, and all interactions of these 
discontinuities are evaluated by locally exact solu- 
tions. An example of the code applied to the flow 
field produced by a series of diamond shapes inside a 
duct (simulating a scramjet combustor) is shown in 
figure 12 (from ref. 13). Included in the figure are the 
shock waves (heavy lines), the vortex sheets (dashed 
lines), and the isobar pattern. Some discontinuities 
have been eliminated, because they were found to be 
extremely weak. The effort throughout this work is 
to provide the designer and the analyst with an 
accurate, reliable, user-oriented tool. 

Exhaust 

The afterbody, like the forebody, of an aircraft 
with an integrated propulsion system must be de- 
signed in light of both engine thrust and vehicle 
aerodynamics. The engine exhaust flow, because of 
physical area limitations, may be underexpanded at 
the nozzle exit, and, to obtain maximum propulsive 
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Figure 12. -Simulated seramjet combustor (from Salas, ref. 13). 


efficiency, the vehicle afterbody undersurface is used 
to provide additional expansion. This results in a 
three-dimensional nozzle flow the boundaries of 
which are defined both by the solid boundary of the 
nozzle wall and by the boundary separating the 
nozzle flow from the vehicle external flow. The. 
dominant features of the seramjet nozzle are as 
follows: 

(1) The gas mixture is composed of hydrogen/air 
combustion products. 

(2) The flow at the combustor exit has nonuni- 
form composition, variable stagnation properties, and 
a complex wave field. 

(3) There is a rapidly expanding internal nozzle 
flow field with waves generated by and reflected off 
multiple surfaces. 

(4) There is interaction of the nozzle exhaust flow 
field with the nonuniform vehicle external flow as 
well as intermodule interactions. 

Such flows for supersonic aircraft (see fig. 13, ref. 14) 
are being studied by Advanced Technology Labora- 
tories, Inc., and NASA Lewis Research Center. 

Two programs developed for flow field calcula- 
tions are CHAR3D and BIGMAC. CHAR3D is a 
reference plane characteristic code with wave- 
preserving network, nonisentropic pressure density 
relation, and conservationlike cross derivatives. 
BIGMAC is a reference plane finite difference code 
employing conservation variables and a one-sided 



Figure 13. -Exhaust How field for the hypersonic research 
airplane (from Dash and Del Guidice, ref. 14). 

difference algorithm. Common features of both pro- 
grams include the following: 

(1) Reference plane grid network with respect to 
three-coordinate systems (Cartesian, cylindrical, and 
line source) 

j (2) Quasi-streamlines followed in reference planes 

(3) Reference characteristic calculation at all 
boundary points 

(4) Dual cubic spline geometry package 

(5) Equilibrium chemistry package 

(6) Discrete treatment of plume boundary and 
underexpansion shock 

An example of a calculation for a three-dimensional 
nozzle exhaust illustrated in figure 14 is shown in 
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Figure 14.-F-15 with proposed rectangular nozzle. 



x/d 


Figure 15. -A simple model of exhaust from a three-dimensional rectangular nozzle. 


figure 15. This calculation is the result of an ADI 
technique applied to the parabolic Navier-Stokes 
equations. 

CONFIGURATIONS 

The following discussion on configurations is given 
in terms of two-dimensional airfoils, subsonic and 
supersonic aircraft, transonic aircraft, and Space 
Shuttle applications. 

Airfoils 

Current computer programs for the computation 
of two-dimensional, steady, in viscid aerodynamic 


flows over single element airfoils are fast and flexible, 
having been highly developed over the past few years. 
As such they provide a valuable analysis tool for 
subsonic, transonic, and supersonic flow. A summary 
of this research is given in table 5. 

The approach most widely used to account for 
viscous effects in an analysis was first suggested by 
Prandtl (ref. 15). The boundary layer displacement 
thickness is added to the original geometry and 
produces an equivalent inviscid shape that represents 
the displacement of the inviscid flow streamlines by 
the boundary layer. A number of researchers have 
coupled an inviscid code with a boundary layer code 
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TABLE 5—2-Dimensional Configurations 


Investigator 

Method 

Application 

Code status 

R. L. Barger and C. W. 
Brooks, Jr. 

Finite difference, inviseid 

Airfoil design and analysis 

Operational; plans include 
extension to 3 dimen- 
sions, rotor blades, and 
engine inlets 

H. L, Morgan 

Integral method: coupled 
boundary layer/potential 
flow 

Subsonic airfoil design and 
analysis 

Operational for analysis of 
multielement airfoils 

E. Murman 

Finite difference, inviseid 

Subsonic and transonic air- 
foil design and analysis 

Operational; plans include 
extension to full speed 
range and coupling with 
boundary layer code 

R. F. Warming and R. M. 
Beam 

Finite difference, hybrid 
scheme, inviseid 

Dynamic airfoil analysis 

Operational 

E. D. Martin 

Semidirect finite difference, 
inviseid 

Transonic thin airfoil analysis 

Operational; plans include 
extension to arbitrary air- 
foil 

R. Hicks 

Numerical optimization, 
inviseid 

Subsonic and transonic air- 
foil design 

Operational; plans include 
developing user package 
and extension to 3- 
dimensional wings 

G. S. Deiwert 

Explicit finite difference, 
time-marching Navier- 
Stokes 

Transonic airfoil analysis 

Operational 

P. Bavitz 

Coupled boundary layer/ 
potential flow 

Transonic airfoil analysis 

Operational 

L. Olson 

Coupled boundary layer/ 
potential How 

Multielement airfoil tran- 
sonic analysis 

Operational; plans include ex- 
tension to 3 dimensions 

F. Bauer, P. Garabedian, and 

Complex characteristics in 

Transonic shockless airfoil 

Operational; plans include 

D. Korn 

hodograph plane 

analysis 

developing design 
capability 

A. Jameson 

Finite difference, relaxation, 
with fast elliptic solver 

Transonic airfoil analysis 

Operational 

L. A. Carlson 

Finite difference, inverse 

Transonic airfoil design and 
analysis 

Operational 

J. E. Carter and S. F. 
Womom 

Coupled boundary layer/ 
potential flow 

Calculations through small 
separation regions in air- 
foil analysis 

Operational 

T. J. Coakley and J. G. 
Marvin 

Explicit finite difference, 
Navier-Stokes 

Shock-boundary layer inter- 
action 

Operational 

C. M. Hung 

Explicit finite difference, 
Navier-Stokes 

Compression corner 

Operational 

S. H. Goradia 

Integral method: coupled 
boundary layer/potential 
flow 

Subsonic airfoil 

Developmental 

U. B. Mehta 

Implicit finite difference, 
Navier-Stokes 

Unsteady separation 

Operational 

J. F. Thompson 

Finite difference, Navier- 
Stokes 

Incompressible multielement 
airfoil 

Operational; plans include 
adding compressibility 
effects 

J.C. Wu 

Integro-differential, Navier- 
Stokes 

Unsteady separation on air- 
foil 

Operational 


to provide a detailed analysis of the interacting flow 
field. (These codes have also been used successfully 
for shock/boundary layer investigations providing the 


interaction does not cause a large separation bubble 
to occur.) Results given in figure 16 (from ref. 16) 
show the complete adequacy of the analysis method 
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Figure 16.-GA(W)-1 airfoil comparison (from Barger and Brooks, ref, 16). (c m = mass coefficient.) 


for subsonic airfoils at angle of attack below stall. At 
higher angles of attack, where there is flow separa- 
tion, the theory breaks down. Viscous theory is 
sometimes inadequate before separation occurs. Thus, 
as illustrated in figure 17 (from ref. 17), the theory 
cannot be depended upon to predict the characteris- 
tics of very thick airfoils (e.g., 21 percent) at low 
Reynolds numbers where large viscous effects domi- 
nate. 

For analysis of multielement airfoils with high-lift 
capability, present analysis programs are not as 
accurate as for single-element airfoils although they 
are useful engineering analysis and design tools for 
high-lift devices. The two-element airfoil configura- 
tion shown in figure 18 (from ref. 18) and the 
four-element airfoil configuration shown in figures 19 
(from ref. 19) and 20 (from ref. 20) are illustrative of 
the analysis capability. 


Current transonic analysis methods such as those 
developed at New York University (by Garabedian 
(ref. 21) and Jameson (ref. 22)) and by Murman (ref. 
23) effectively handle the analysis of transonic 
airfoils. This analysis capability is shown in figure 21 
for the case of a supercritical airfoil and in figure 22 
for the case of a conventional airfoil (ref. 24). Of 
particular interest is the good prediction of shock 
location and strength. 

Numerical methods for airfoil design are also 
highly developed for subsonic and supersonic flow 
and are in the process of development for transonic 
flow. The best-known and most widely used design 
methods are the hodograph methods of Bauer, 
Garabedian, and Korn (ref. 25). While these methods 
yield accurate results, they involve complex mappings 
and transformations. Another design technique is to 
use a direct method to analyze the flow about a given 
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Figure 17. -Analysis of general aviation airfoil (from McGhee and Beasley, ref. 17). M - 0.15;/? = 3.0 

X 10 6 ; fixed transition. 
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Figure 18.- Analysis of two-element airfoil with single slotted flap (from Olson and Dvorak, ref. 18). 
5 flap = 30° * fa) Diagram. ( b ) Section pressure distribution, a = 7.5°. (c) Section lift. 
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Figure 18 (concluded).-(tf) Lift 



Flap gap (percent chord) 
-O Experimental 

(e) -n-Theoretical 

is drag. ( e ) Flap gap optimization. 



0 “ 5°/ 1 5° 

□ 147-33° 

a 30765° (flow is separated 
on the aft flap element) 

Theoretical 

Figure 19. -Analysis of four-element airfoil (from Morgan, 
ref. 19). 5 V = 6 vane ; 5/ = 6 flap ; = -45°. 



Figure 20. -Analysis of four-element airfoil (from Goradia 
and Lilley, ref. 20). 6, = -45° ; 6 V = 14° ; b f = 33°. 


airfoil and then, based upon this result, to modify the and a large number of iterations, but it can give good 

airfoil to satisfy the design conditions. This usually results. Research is currently underway at NASA 

requires extensive experience on the part of the user Ames Research Center (Hicks, ref. 26) in which a 
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Theory; c f = 0.719 

o Experiment; c n =0.712 


Theory; c t - 0.248 

o Experiment; c n = 0.242 


Figure 21. -Analysis of Korn supercritical airfoil (from 
Bavitz, ref. 24 ). M = 0.752 ; R = 20.95 X 10 6 ;a = 1.25°. 


Figure 22. -Analysis of conventional airfoil NACA 65 j-2 1 3 
(from Bavitz, ref. 24). a = 0.5 ; 4/ = 0.753; R = 52.60 X 
10 6 ; a = -0.10°. 


small perturbation analysis program is combined with 
an optimization program based on the method of 
feasible directions. This design capability extends to 
the optimization of various characteristics, as illus- 
trated in figure 23 (from ref. 27) for supersonic flow 
and in figure 24 (from ref. 28) for transonic flow. A 
third design method is an inverse method in which 
the airfoil surface pressures or velocities are specified 
and the airfoil shape is determined from the calcula- 
tion. This approach requires a knowledge of what 
would be a good pressure distribution. While the user 
can select a pressure distribution that would result in 
a given lift and moment and satisfy other flow 
characteristics, the resultant airfoil design may or 
may not be physically and/or structurally reasonable. 
Current research using the inverse method is con- 
cerned with providing improved airfoils of practical 
design. 

The power of the theory for the analysis and 
design of airfoils, as opposed to past empirical 
methods that involved scaling of ordinates to create 


airfoil families, is evident. However, current theories 
have been shown to be inadequate for use in the 
analysis and design of airfoils in situations wherein 
significant amounts of flow separation occur. Until 
reliable turbulence modeling procedures or far more 
powerful computers are developed, it will be neces- 
sary to obtain stalled airfoil characteristics by means 
of careful experimentation. However, several analyti- 
cal investigations underway show promise of being 
able to at least qualitatively predict the features of 
two- and three-dimensional bodies with extensive 
flow separation. A Navier-Stokes calculation of an 
airfoil undergoing unsteady laminar separation is 
shown in figure 25 (from refs. 29 and 30) to 
reproduce the major flow features revealed in wind 
tunnel flow visualization studies. 

A different approach to flow separation employs a 
wake description of discrete vortexes arising from the 
separation of shear layers at the surface. The discrete 
vortexes connect and diffuse downstream to form an 
unsteady wake. Figures 26 (from ref. 31) and 27 
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two-dimensional shapes at angle of attack, and 
programing is underway to extend the method to 
arbitrary three-dimensional bodies. 

Subsonic and Supersonic Aircraft 


(a) 



(b) 


Figure 25. -Leading edge stall, (a) Angle of attack = 15°; 
Re = 10 3 ; 9 percent thick symmetrical airfoil; from 
Mehta and Lavan (ref. 29). ( b ) Wind tunnel flow visuali- 
zation (from Prandtl, ref. 30). 


Fluctuating 

Separation line vortical flow 



Boundary layer 


Figure 26. -Schematic drawing of three-dimensional un- 
steady, separated flow (from Leonard, ref. 31). 



(from ref. 32) show the method as applied to a sphere 
and indicate that the calculations produce vortex 
flow patterns similar to those observed experimen- 
tally. Analyses have been completed for a variety of 


At subsonic and moderate supersonic velocities 
and lower altitudes the airflow around a configura- 
tion remains in chemical equilibrium, and design 
information can be obtained from the perfect gas 
flow calculations and wind tunnel testing of appropri- 
ately scaled models. Current three-dimensional ana- 
lytical methods for complete configurations are usu- 
ally based on potential methods and are suitable only 
for inviscid flow. (See table 6.) The highly developed 
Boeing program has made extensive use of the 
paneling procedure for a variety of configurations. 
For the aerodynamic representation, the configura- 
tion surface is divided into a number of panels. A 
constant source distribution is then imposed on the 
body panels and a linear vortex distribution is used to 
represent the wing and tail panels. Typical results for 
wing pressures for two types of wing/body configura- 
tions (a supersonic transport and a maneuvering 
fighter) are presented in figure 28 (from ref. 33). 
Current research is directed toward incorporating the 
subsonic doublet-lattice and vortex spline aero- 
dynamics technology to improve programing accu- 
racy and efficiency. 

For high-speed aircraft, the effects of separation- 
induced vortex flows are becoming increasingly im- 
portant, particularly optimization of vortex-induced 
maneuver lift. Prediction of lift and pressure distribu- 
tions for delta wings at angle of attack have been 
attempted with linear lifting surface theory, slender 
wing theory, and the Smith conical method. How- 
ever, these approaches have been inadequate to 
handle the vortex flow emanating from the leading 
edge of such wings. Polhamus’ recent leading-edge 
suction analogy method (ref. 34), which has proven 
to be extremely successful in predicting the lift on a 
number of configurations in which vortex tlow is a 
significant contributor, does not predict pressure 
distributions. A complete three-dimensional paneling 
solution, which has recently become available from 
Boeing and which accounts for the free-sheet vortex 
and frozen wake, accurately predicts such pressures. 
This free-sheet method is compared with data and 
some earlier theories in figure 29 (from ref. 35). 

Some application of the viscous equations for 
steady three-dimensional flow is now being made 
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TABLE 6 -Subsonic and Supersonic Aircraft 


Investigator 

Method 

Application 

Code status 

L. Morino 

Finite element 

Subsonic and supersonic con- 
figuration design 

Operational; plans include 
extension to transonic 
How 

W. D. Middleton 

Linear theory 

Integrated system design and 
optimization 

Operational for supersonic 
flow' code currently being 
optimized 

P. Kutler and L. Sakcll 

Finite difference, invisdd 

Shock-on-shock interaction 

Operational 

J. A. Weber, G. W. Brune, 
and F. T. Johnson 

Paneling 

Wing leading edge 

Operational for subsonic 
How 

J. E. Lamar 

■ Suction analogy 

Augmented vortex lift inter- 
action effects 

Operational; plans include 
extension to round lead- 
ing edge wings 

P. E. Rubbert, F. T. 

Johnson, and F. E. Ehlers 

i Paneling and doublet spline 
fitting 

Full configuration design 

Design operational; plans in- 
clude extensions for 
dynamic loads and active 
control configurations 

A. Leonard 

Coupled boundary layer/ 
potential flow 

Interacting vortex filaments 

Operational 

E. Lan 

Wing plus thick jet theory 

Jet exhaust/wing interaction 

Operational for subsonic- 
flow with arbitrary jet 
shape and deflection; 
plans include extension 
to higher much numbers 

V. Rossow 

Lagrangian marching 

Vortex alleviation in incom- 
pressible flow 

Operational 

M. L. Lopez and C. C. Shcn 

Jet flap theory and elemen- 
tary vortex distribution 

Propulsive lift on wings 

Operational; plans include 
extending calculation to 
full body 

J. R. Tulinius, D. S. Miller, 
J. L. Thomas, and R. J. 
Margason 

Paneling plus jet Hap 

Full configuration 

Developmental; will be used 
for configuration design 

S. C. Lubard 

Fully implicit finite dif- 
ference, parabolic Navier- 
Stokes 

Analytic configuration in 
supersonic flow 

Operational 

S. G. Rubin and C. L. Lin 

Two-step finite difference, 
parabolic Navier-Stokes 

Analytic configuration in 
supersonic flow 

Operational 

R.S. Hirsh 

Finite difference (ADI) 

3-dimensional nozzle ex- 
haust 

Operational for supersonic 
flow; plans include ex- 
tension to subsonic 

W. S. Llewellen and R. D. 
Sullivan 

Finite difference, x- 
m arching viscous 

Subsonic wake flow (now 
being applied to vortex 
mixing) 

Operational 


where diffusion in -one direction can be neglected. 
The governing equations, the so-called parabolic 
Navier-Stokes equations, have been applied to ana- 
lytic configurations in supersonic flow, supersonic 
three-dimensional nozzle exhausts, and subsonic wake 
flow, and are currently being applied to subsonic 
vortex mixing. 

The development of computational methods capa- 
ble of accounting for various jet interaction effects is 
of interest because of the need to analyze highly 


integrated propulsion system/airframe concepts. This 
is a particularly vital area in vertical takeoff and 
landing (VTOL) research and in the upper-surface 
blowing (USB) propulsive-lift concept, which is prom- 
ising for high-performance aircraft with an inherently 
lower noise level because the engine nozzle and 
exhaust are located above the wing and thereby 
shielded by the wing with respect to noises radiated 
to the ground. Research in this area involves applica- 
tion of wing plus thick jet theory and elementary 
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Constant source distribution 
on body panels 



(a) 





O Experiment, 
upper 
surface 

□ Experiment, 
lower 
surface 
— Theory, 
planar 
boundary 
condition 


Figure 28. -Pressure calculations from paneling method of 
Ehlers, Johnson, and Rubbert (ref. 33). (o) Aerody- 
namic representation. ( b ) Supersonic transport wing. M 
= 0.60; a = 5.63°. (c) Supersonic transport wing. M = 
1.90; a = 5.28°. ( d ) Fighter wing; large camber. M = 
0.40; a - 8.37°. (e) Fighter wing; moderate camber. M 
= 0.40; a = 3.93°. (/) Fighter wing; no camber. M = 
0.26; a = 6.4°. 



Theories 

o,o,a,0 Experiment — Attached flow 

Conical vortex flow 

Free sheet method 

Figure 29. -Edge separation vortex flow (from Weber et al., 
ref. 35). A = 1.46; a = 14°;Af - 0. 


vortex distribution. An example of the predictive 
capability of the thick jet theory in estimating the lift 
and pitching moment on a USB configuration is 
shown in figure 30 (from Lan, ref. 36). It can be seen 
that this theory is a distinct improvement over the 
thin jet flap approach that neglects the wing/jet 
interaction effect associated with a thick jet. Dr. Lan 
is currently extending the theory to account for a 
round jet blowing over the wing, to examine the 
relative importance of jet entrainment and inviscid 
wing/jet interaction effects for configurations at high 
subsonic speeds, and to account for fighter aircraft 
designs that have integrated propulsion systems. 
Unlike the USB concept, a number of propulsive lift 
approaches involve the use of jet flaps that do not 
interact with the wing. The elementary vortex distri- 
bution (EVD) jet flap theory is shown to be 
particularly effective in predicting the characteristics 
of such configurations. An example of such capability 
is shown in figure 31 (from ref. 38). 

Multiple trailing vortices from large aircraft con- 
tinue to pose a serious hazard for following aircraft. 
To provide further insight into the rollup process and 
the mechanisms for eventual breakup of wake vor- 
texes, to verify existing wake vortex predictive 
techniques, and to supplement experiments, un- 
steady, three-dimensional numerical simulations of 
the flow field using vortex filaments are being 
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I 


] 0.1667 




o Experiments {ref. 37) 

— Present thick jet theory 
-- Present thin jet flap theory 


Figure 30. -Jet flap theory for USB configuration (from Lan, 
ret. 36). (Sketch of USB configuration appears in ref. 37.) 
The center of gravity is 0.49 below the wing plane. 


performed. The multiple vortexes can merge to form 
a pair of concentrated vortexes with high velocities or 
can reduce to two diffuse vortexes with low veloc- 
ities. A study of vortex merging is underway to 
investigate this phenomenon and analyze the possibil- 
ity of controlling it. One investigation at NASA Ames 
(ref. 31) uses a Lagrangian marching procedure. 
Another simulates the wake vortex flow field by a 
numerical method based on an Eulerian (fixed grid) 
finite difference approach (fig. 32). Both the incom- 
pressible inviscid and viscous and the compressible 
inviscid equations for unsteady, two-dimensional flow 
are being studied as to the relative merits of each for 
determining the wake flow field. The algorithms for 
both the compressible and incompressible equations 
are second-order accurate in time and second- or 
fourth-order accurate in space. The method used for 
the incompressible equations is semi-implicit, while 
the method used for the compressible equations is 
explicit. After a decision has been made as to which 
formulation is best, that procedure will be extended 
to include all three space dimensions. Results of 
integrating the compressible Euler equations 0.73 
span downstream in a wake composed of four 
vortexes are shown in figure 32(c). 

Another area where numerical simulations are of 
interest is the prediction of flow characteristics over 
missiles and missile-shaped configurations. A method 
has been developed for solving the exact nonlinear 
potential equation for axisymmetric flow about blunt 
or pointed bodies with surface anomalies such as 
sharp-cornered boattails and flares. The method 
allows consideration of quite general shapes and is 
operational throughout the speed range. In the results 
shown in figure 33 (from ref. 39) for a subsonic 
freestream, the flow becomes supercritical in the 
shoulder region, shocks down on the cylindrical 
portion, expands again to supersonic at the sharp 
boattail convex corner, and is shocked down again. 
The calculated flow compresses again to nearly 
stagnation pressure in the concave corner, which is 
obviously a poor simulation of the wake separation. 
Figure 34 (from ref. 39) shows results of a calculation 
for a typical launch vehicle shape in supersonic flow 
with a simulated exhaust plume. 

Transonic Aircraft 

Most modern high-speed aircraft achieve optimum 
cruise performance and maneuverability when flying 
in the transonic regime where the flow is extremely 





30 


A SURVEY OF COMPUTATIONAL AERODYNAMICS IN THE UNITED STATES 



(b) (c) T= 0.73 spans 

Figure 32. -Wake vortex problem (from Paul Kutler, unpublished), (a) Diagram. (6) Accuracy check 
for incompressible, inviscid vortex, (c) Multiple vortex pattern (compressible equations). 



Figure 33. -Pressure distribution for ogive/cylinder/boattail Figure 34. -Pressure distribution for staged vehicle with 
(from Keller and South, ref. 39). M = 0.83; 193 X 49 simulated plume (from Keller and South, ref. 39). M - 

mesh. 1.2; 193 X 49 mesh. 
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TABLE 1 .—Transonic Aircraft 


Investigator 

Method 

Application 

Code status 

P. A. Newman and E. B. 
Klunker 

Finite difference, inviscid, 
small disturbance 
approximation 

Full configuration 

Operational 

F.R. Bailey and W. F. 
Ballhaus 

Finite difference, relaxation, 
inviscid, small disturbance 
approximation 

Wing/body interaction 

Operational 

J. C. South, Jr., and J. D. 
Keller 

Finite difference, relaxation, 
mapped coordinates, full 
potential 

Missile shapes 

Operational 

R. W. Barnwell 

Finite difference, relaxation 

Wing/body interaction 

Operational 

P. R. Garabedian and D. 
Korn 

Complex characteristics, full 
potential 

Wings 

Operational 

A. Jameson 

Finite difference, relaxation | 
1 

Swept and yawed wing 

Operational 


sensitive to aircraft configuration. Because of the 
special problems associated with both cruising and 
maneuvering flight at near sonic speed, a significant 
effort is being made to develop theoretical transonic 
flow analysis methods. (See table 7.) The highly 
nonlinear nature of the equations makes linear 
analysis of little value; consequently, solutions to the 
nonlinear governing equations must usually be ob- 
tained. 

The status of the present effort aimed at com- 
puting transonic flows over wing/body configurations 
by solving finite difference approximations to the 
transonic small-disturbance equation is illustrated in 
figures 35 to 39 (from refs. 41 and 42). The 
conservative finite-difference equation is derived by 
applying the divergence theorem to the conservation 
integral of the governing equation over a computa- 
tional cell. Such an approach assures that the shock 



Figure 35. -Comparison of flight, wind tunnel, and computed 
pressure distribution on upper surface of C-141 wing from 
Lomax, Bailey, and Ballhaus (ref. 41). = 0.825. 


jump relations implied by the governing equation are 
obtained in the finite-difference approximation. Com- 
parisons of computed and experimental pressure 
distributions for the C-141 wing shown in figure 35 
illustrate the accuracy of the technique in predicting 
the shock strength and location. (The comparison 
also shows that, in this case at least, inviscid 
calculations predicted the flight situation more accu- 
rately than low Reynolds number wind tunnel tests.) 
Results are shown in figures 36 through 38 for several 
isolated wings and in figure 39 for a wing/body 
configuration. In these figures, the two types of 
computed results are fully conservative relaxation 
(FCR), which models the true inviscid shock jump, 



Experiment: RAE, Re. - 1.4 X 10 6 , 
8- by 6 ft tunnel 
Computed: NCR and FCR 


Figure 36. -Pressure distribution for Royal Aircraft Establish- 
ment (RAE) wing C (from Bailey and Ballhaus, ref. 42). 
Afoo = 0.95; a = 0° ; sweep angle A C f 2 = 44°; aspect ratio 
(AR) = 3 6; taper ratio (TR) = 0.333; RAE 101 section; 
5.4 percent streamwise. (FCR = fully conservative relaxa- 
tion; NCR = nonconservative relaxation; Re c = Reynolds 
number based on chord.) 
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□ Lower surface 
Theory, small disturbance 
_ FCR 
__ NCR 

Figure 37. -Pressure distribution for the ONERA M6 wing (from Bailey and Batlhaus, ref. 42). = 

0.84; a = 3.0°; A LF = 30°; AR = 3.86; TR = 5.83; ONERA D section; 9.8 percent streamwise. 


and nonc'onservative relaxation (NCR), which obtains 
a modified shock jump that often agrees better with 
experimental results because of viscous effects. 

Tunnel wall ; merference plays an important part 
in transonic flow experiments; therefore, inclusion of 
tunnel wall boundary conditions is necessary in all 
transonic calculation codes before valid comparisons 
with experimental data can be made. Calculations to 
simulate a number of conventional tunnel wall 
boundary conditions have recently been made using a 
relaxation technique to solve a small disturbance 
potential equation. Initial numerical results for finite 
lifting wings show that the embedded supersonic 


bubble, generally terminated by a shock wave, is 
distorted compared with that for free air. Computed 
tunnel wall pressure distributions for several con- 
ventional linear boundary conditions compared with 
the free air results at the wall location are shown in 
figure 40 (from ref. 43). The present nonlinear finite 
difference formulation of the problem affords far 
more flexibility with respect to the tunnel wall 
boundary condition than was possible in conventional 
linear correction theory. 

For some applications, however, approximate 
methods are satisfactory and substantial reductions in 
computer time and storage can be realized. For these 
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Experiment, ONERA 
O Upper surface 
□ Lower surface 
Theory, small disturbance 

FCR 

NCR 

Figure 38. -Pressure distribution for the ONERA M6 wing (from Bailey and Ballhaus, ref. 42). 
0.92; a = 3.0°; A LE = 30° ; AR = 3.86; TR = 0.538; ONERA D section; 9.8 percent streamwise. 


applications a rapid approximate method for calcu- 
lating transonic flow about lifting configurations was 
developed through an analytical study of the effect of 
lift on the transonic area rule. This method reduces 
the three-dimensional fluid mechanics problem to a 
two-variable computational problem. A highly effi- 
cient finite difference relaxation procedure is then 
applied. Calculations using this method to illustrate 
the effect of including tunnel wall boundary condi- 
tions are shown in figure 41 (from ref. 44). This wall 
interference effect is observed as a shift in the shock 
wave and sonic line locations. 

The examples presented in this section on tran- 


sonic aircraft have been computed solutions to the 
small disturbance equation, except for the approxi- 
mate method based on the lifting area rule. Some 
solutions to the full-potential equation (exact com- 
pressible inviscid equation) now exist for simple 
geometric configurations such as missile shapes and 
moderately swept wings. An example of the analysis 
capability of solving the full potential equation for a 
relaxation technique for transonic flow calculation 
was presented in figures 33 and 34. Extending 
transonic flow calculations to include viscous effects, 
either through coupling an inviscid procedure with a 
boundary layer procedure or by solving the Navier- 
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Experiment: 

Langley 8 -ft wind tunnel 
Re, = 2.0 X 10 6 

o Upper surface transition natural 
□ Upper surface transition fixed 
& Lower surface transition natural 
O Lower surface transition fixed 
NCR calculation 


y/b = .95 
y/b = .80 

y/b = .60 



y/b = .40 


Figure 39. -Pressure distribution for swept wing/fuselage configuration (from Bailey and Bailhaus, ref. 
42). = 0.93; a = 0° ; A C/4 = 45° ; AR = 4; TR = 6; NACA 65 A006 streamwise. 



Free air 

Solid wall 

Porous wall 

Slotted (S = 0.104) 

Slotted upper and 

lower solid side 
(S = 0.104) 

Contoured upper and 

lower straight solid 
side 

Figure 40.— Tunnel wall pressure distributions in the wing 
root plane (from Newman and Klunkcr, ref. 43). NACA 
63 A006 section ; rectangular plan ; AR = 32 /9 ; M = 0.9 3 ; a 
= 2.76° ;/i/c = 5.3. 

Stokes equations for transonic flow over configura- 
tions, is not feasible without a more detailed under- 
standing of the transonic flow regime. Both detailed 


experimental data and numerical results including 
tunnel wall effects are required to determine a 
satisfactory means of approximating the effects of 
complex boundaries. 

Space Shuttle 

At high velocities and altitudes, significant non- 
equilibrium effects will exist over a large area of the 
Space Shuttle vehicle. The calculation of flow with 
finite rate chemistry is, therefore, needed by the 
Space Shuttle designer as a source of information 
relative to the prediction of heat transfer rates, 
boundary layer effects, and aerodynamic loads acting 
on the aircraft. Complex chemically reacting flow 
phenomena cannot be scaled, and because current 
test facilities are incapable of performing full-scale 
experiments in this regime, the designer must turn to 
the computer solution for a realistic description of 
the flow field at actual flight conditions. The major- 
ity of the present calculations for the Space Shuttle 
involve the solution of the inviscid equations. For 
previous manned atmospheric entry vehicles, non- 
equilibrium effects could be determined by boundary 
layer analysis alone, using equilibrium properties for 
the inviscid external flow. For the Space Shuttle, 
however, a large part of the inviscid flow will also be 
nonequilibrium; therefore, a fully reacting inviscid 
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Free air flow Flow in tunnel 



(b) 



— Reacting gas 
— Perfect gas 

(a) 



— Reacting gas 
_ . Perfect gas 
(b) 


Figure 41. -Tunnel wall effects in transonic flow (from 
Barnwell, ref. 44). {a) Wing plane results for free air and 
slotted tunnel. M = 0.98; a = 12°. {b) Effect of tunnel 
wall on shock wave and sonic line locations. M ^ - 1.02; a 
= 10 °. 

solution is needed to provide the necessary boundary 
conditions for a solution of the reacting boundary 
layer. 

Blunt body flows develop a thin entropy layer 
near the body surface in which strong gradients of 
flow properties occur. This entropy layer occurs also 
in reacting flows but is more complex because of the 
chemical species gradients. Bow shock shapes calcu- 
lated using a code based on the method of character- 
istics are shown for the side and plan views and for 
two cross sections in figure 42 (from ref. 45). To 
provide some feel for the effects of nonequilibrium 
chemistry, the Space Shuttle results computed incor- 
porating real gas effects are compared with perfect 
gas (frozen) computations at the same free stream 



— Reacting gas 
Perfect gas 

(c) (d) 


Figure 42. -Shock shape for Space Shuttle 147 body (from 
Rakich and Kutler, ref. 45). a = 30°; V ^ = 6.7 km/s; 
altitude = 65.5 km. (j) Side view. ( b ) Plan view, (c) 
Cross section at x = 4.75 m. ( d ) Cross section at x = 
15.2 m. 
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mach number. The shock standoff distance is de- 
creased by the chemical reactions. Note that the 
location of the point of bow shock impingement on 
the wing is important to the heat shield design. 

At moderate angles of attack (up to approximately 
30°), the subsonic portion of the flow over winged 
lifting vehicles, such as the Space Shuttle, is confined 
to the nose region. An effort at Grumman has 
resulted in the development of an inviscid computer 
code for solution of the moderate angle of attack 
problem (30° to 35°) that provides inputs for a 
three-dimensional boundary layer code. All shock 
waves are computed explicitly as discontinuities in 
the flow. Figure 43 (from ref. 46) shows a compari- 
son of calculated and experimental pressure distribu- 
tions along the windward and leeward symmetry 
plane of a Space Shuttle configuration at angle of 
attack. The agreement is generally good even on the 
leeward side, except in the nose region where the 
difference is thought to be due to an oversimplifica- 
tion of the vehicle geometry for the calculations. 
More accurate geometric modeling will solve this 
problem. 

Calculations for the Space Shuttle configuration 
using both a three-dimensional inviscid code coupled 



* = 55° 



— Theory 

o Experiment at Ames 3.5 ft 
hypersonic wind tunnel 

Figure 4 3. -Longitudinal surface pressures in the Space 
Shuttle body (from Marconi, Yaegcr, and Hamilton, ref. 
46). M - 7.3 ; 7 = 1.4;a = 30 °\L = 107.8 ft. 


with a boundary layer code and a three-dimensional 
Navier-Stokes code were made. Computed and experi- 
mental heating rates are compared in figure 44 (from 
ref. 47) for the windward and leeward centerlines and 
for one transverse section. As illustrated, the three- 
dimensional Navier-Stokes solution is extremely “user 
sensitive” in regard to finite-difference grid selection, 




Figure 44. -Heating rate on Space Shuttle (from Goodrich et 
al., ref. 47). ( a ) Windward centerline. ( b ) Leeward 
centerline, a ~ 30°. (See key on next page.) 
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Computed data: 

— Newtonian/boundary layer (normal shock entropy) 
Newtonian/boundary layer (parallel shock entropy) 
UK Three-dimensional Navier-Stokes 
(grid influence) 


Unpublished wind tunnel data: 

o AEDC - B, = 7.9 ) 

□ AEDC - F, = 12 > Re^ L * 10 6 

§ Calspan, 96-in.; M x = 7 to 12/ 


(0 


Figure 44 (concluded), -(c) Transverse section, x/L = 0.1. 


TABLE S. Space Shuttle 


Investigator 

Method 

Application 

Code status 

A. W. Rizzi 

Time-split finite difference, 
inviscid 

Nonequilibrium flow, 
chemically reacting 

Operational 

J. V. Rakich and E. B. Pegot 

Coupled boundary layer, 
inviscid 

Windward side heating pre- 
diction, laminar flow, 
perfect gas 

Operational 

C. P. Li, C. K. Houston, and 
R. M. Meyers 

Coupled boundary layer/ 
potential flow, finite dif- 
ference, Navier-Stokes 

Full configuration with real- 
gas effects 

Operational 

F. Marconi and L. Yaeger 

Finite difference, inviscid 

Full configuration with real- 
gas effects 

Operational 

P. Kutler 

Finite difference, inviscid 

Unsteady shock interaction 

Operational 


a feature that discourages its use as a design tool at 
the present time. There is currently a need for more 
physical and numerical analyses in addition to more 


computer power to aid in resolving the difficulties of 
three-dimensional Navier-Stokes calculations. (See 
table 8 for a summary of this research in this area.) 


CONCLUDING REMARKS 


As illustrated in table 9 and figures 45 and 46, the 
cost of numerical simulation continues to decrease as 
computer speed and storage continue to increase; 
consequently, simulations of practical problems of 
interest to the aerospace community are becoming a 
reality. Computer codes to calculate solutions of the 
inviscid approximations and turbulent boundary layer 
approximations to the Navier-Stokes equations are 
highly developed and are currently being applied in 
aircraft design. Current research is directed toward 
developing numerical procedures, turbulence trans- 
port models, and computer codes for obtaining 
solutions to the time-averaged Navier-Stokes equa- 
tions. 


TABLE 9 —Arithmetic Instruction Timings in 
Seconds 


Computer 

Operation 

+/- 

X 

I 

BABBAGE (1800) 

1 

60 

60 

MARK 1 (1937) 

.3 

6 

11.4 

ENIAC (1946) 

2 X 10" 4 

3 x 10 -3 

6 x 10" 3 

CDC 6600 (1963) 

4 X 10- 7 

5 X 10- 7 

2 X 10“‘ 

ILLIAC IV 

6 X 10" 5 

7.8 X 10- 9 

4.9 X 10 _ * 

STAR 100 

2 X 10“' 

4 X 10-' 

4 X 10-' 


Estimated minimum requirement for 


practical computation of turbulent 
eddies in three-dimensional aerodynamic 



Year new computer available 
Figure 46. -Trends of computer speed. 

Numerical solutions of the time-dependent 
Navier-Stokes equations for a mesh point distribution 
sufficiently refined to represent turbulent motion are 
not feasible with current and projected computer 
systems (storage and processing time requirements); 
however, solutions of these equations with suitable 
subgrid scale turbulence modeling is currently feasible 
and should provide insight into the mechanics of 
turbulent flow. As indicated in figure 46, such 
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Figure 45. -Trend of computation cost for computer simulation of a given flow. 
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solutions would require an increase in computational 
soeed of three orders of magnitude -an increase that 
appears to be feasible in view of current develop- 
ments in computer technology' and numerical proce- 
dures. 

The increasing costs of developing and maintaining 
computer codes for numerical simulation currently 
tend to offset the decreasing cost of computer 
processing; consequently, each research objective 
must be carefully analyzed to assure optimum use of 
available resources, both manpower and computer 
systems, if accurate and efficient computer codes are 
to be developed that will have an impact on practical 
fluid mechanics problems. At the same time, we must 
insure that the right kind and amount of experi- 
mental data are obtained to verify and to provide 
necessary empirical inputs to the numerical simula- 
tion to minimize Professor Bradshaw’s “fact gap” 
(fig. 47). 
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